#include  <fftw3.h>
#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>

int main( int argc, char** argv )
{
	/* Initialize */
	int N = 256;
	double PI = 4.0 * atan(1.0);
	double *f, *u;
	double *xGrid;
	fftw_complex *FTf, *FTu;
	fftw_plan plan_r2c, plan_c2r;
	int j, k;
	double multtmp;

	/* Setup */
	xGrid = (double*)malloc(sizeof(double) * N );
	u = (double*)malloc(sizeof(double) * N );
	f = (double*)malloc(sizeof(double) * N );
	FTf = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1) );
	FTu = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1) );
	plan_r2c = fftw_plan_dft_r2c_1d( N, f, FTf, FFTW_ESTIMATE);
	plan_c2r = fftw_plan_dft_c2r_1d( N, FTu, u, FFTW_ESTIMATE);

	/* Compute */	
	for( j = 0; j < N; j++ )
	{
		xGrid[j] = ( (double)(j) + 0.5 ) / (double)(N);
		f[j] = sin( 2.0 * PI * xGrid[j] );
	}
	fftw_execute(plan_r2c);
	for( k = 0; k < N/2+1; k++ )
	{
		multtmp = FTf[k][0];
		FTf[k][0] = FTf[k][0] * cos( PI * (double)(k) / N ) - 
			FTf[k][1] * sin( PI * (double)(k) / N );
		FTf[k][1] = + multtmp * sin( PI * (double)(k) / N ) + 
			FTf[k][1] * cos( PI * (double)(k) / N );
	}

	FTu[0][0] = 0.0;
	FTu[0][1] = 0.0;
	for( k = 1; k < N/2 + 1; k++ )
	{
		FTu[k][0] = FTf[k][0] / ( 4.0 * PI * PI * k * k );	
		FTu[k][1] = FTf[k][1] / ( 4.0 * PI * PI * k * k );
		multtmp = FTu[k][0];
		FTu[k][0] = FTu[k][0] * cos( PI * (double)(k) / N ) + 
			FTu[k][1] * sin( PI * (double)(k) / N );
		FTu[k][1] = - multtmp * sin( PI * (double)(k) / N ) + 
			FTu[k][1] * cos( PI * (double)(k) / N );
	}
	for ( k = 0; k < N/2+1 ; k++ )
	{
		printf( "(%5.2f, %5.2f)",  FTu[k][0], FTu[k][1] );
	}
	fftw_execute(plan_c2r);
	for ( j = 1; j < N ; j++ )
	{
		u[j] = u[j] / N;
	}
	
	/* Output */
	printf("Error=\n");
	for ( j = 0; j < N ; j++ )
	{
		printf( "%5.2f, ", fabs( sin( 2.0 * PI * xGrid[j] ) / ( 4.0 * PI * PI ) -
														 u[j] ) );
	}
	printf("\n");
	for ( j = 0; j < N ; j++ )
	{
		printf( "%5.2f, ",  sin( 2.0 * PI * xGrid[j] ) / ( 4.0 * PI * PI ) );
	}
	printf("\n");
	for ( j = 0; j < N ; j++ )
	{
		printf( "%5.2f, ",  u[j] );
	}

	/* Free space */
	fftw_destroy_plan(plan_r2c);
	fftw_destroy_plan(plan_c2r);
	fftw_free(FTu);
	fftw_free(FTf);
	free(u);
	free(f);
	free(xGrid);
	return EXIT_SUCCESS;
}

